#include <armadillo>


using namespace arma;

class EP
{
	public:

	EP(void)
	{	
		boost::math::normal_distribution<> *d=new boost::math::normal_distribution<>(0,1); 
		_d=d;
	}
	~EP()
	{
		delete _d;
	}

	void operator()(mat X,double *Yt,double lambda,int npass,double alpha)
	{
		int p=X.n_cols;
		int n=X.n_rows;
		mat mu(p,1);
		mat Q(p,p);
		mat I(p,p);
		I.eye();
	       double *Y=new double[n];	
		for(int i=0;i<n;i++)
		{
			if(Yt[i]==0)
			{
				Y[i]=-1;
			}else{
				Y[i]=1;
			}
		}
		Q.eye();
		Q=lambda*Q;	
		cube Qi(n,p,p);
		Qi.fill(0);
		mat Qmui(n,p);
		Qmui.fill(0);
		mat z(n,1);
		mat ps(n,1);
		mat pmi(n,1);
		mat v(n,1);
		for(int k=0;k<npass;k++)
		{
			for(int i=0;i<n;i++)
			{
				mat oldmu=mu;
				mat oldQ=Q;
				mat temp=Qi(span(i,i),span::all,span::all);
				mat Qmi=Q-temp;
				//mat Lmi=chol(Qmi);
				mat mumi=inv(Qmi)*(Q*mu-Qmui(i,span::all).t());
				mat varz=1+X(i,span::all)*(inv(Qmi)*X(i,span::all).t());
				mat ratz=Y[i]*X(i,span::all)*mumi/sqrt(varz(0,0));
				//cout << "p" << << "\n";
				double pz=cdf(*_d,ratz(0,0));
				mat diffpz=(pdf(*_d,ratz(0,0))/sqrt(varz(0,0)))*Y[i]*X(i,span::all);
				double a=(1/sqrt(varz(0,0)))*(pdf(*_d,ratz(0,0))/pz);
				mat b=pow(a,2)+a*Y[i]*(X(i,span::all)*mumi)/varz(0,0);
				v(i,0)=1/b(0,0)-varz(0,0)+1;
				mat foo=Y[i]*X(i,span::all)*mumi+a/b(0,0);
				pmi(i,0)=foo(0,0);
				mat bar=pz*sqrt(1+(varz(0,0)-1)/v(i,0))*exp(pow(a,2)/(2*b));
				ps(i,0)=bar(0,0);
				mat foo1=inv(Qmi);
				mu=mumi+foo1*(diffpz.t()/pz);
				
				mat diff2pz=(diff_normpdf(ratz(0,0))/varz(0,0))*(X(i,span::all).t()*X(i,span::all));
				mat toto=inv(Qmi)*(diff2pz/pz)+I;
				mat Sigma=toto*inv(Qmi)-(mumi-mu)*((mumi-mu).t());
				Q=inv(Sigma);
				mat T=Qi(span(i,i),span::all,span::all);
				mat ncont=alpha*(Q-Qmi)+(1-alpha)*T;
				Qi(span(i,i),span::all,span::all)=ncont;
		//	cout << Q << "\n";	
				Qmui(i,span::all)=(alpha*(Q*mu-Qmi*mumi)).t()-(1-alpha)*Qmui(i,span::all);
	//		cout << mu.t() <<"\n";		
			}
		
		}
		_m=mu;
		_S=inv(Q);
		//cout << _S;
		delete[] Y;

	}
	double diff_normpdf(double x)
	{
		return -x*pdf(*_d,x);
	}
	mat Get_Sig(void){return _S;}
	mat Get_Mu(void){return _m;}
	private:
	mat _m;
	mat _S;
	boost::math::normal_distribution<> *_d; 
};

